I7l 


RD-A150  926  RESPONSE  OF  SATURATED  SOILS  TO  DYNAMIC  LOADING(U)  OHICT 
STATE  UNI V  RESEARCH  FOUNDATION  COLUMBUS 
R  S  SANDHU  ET  AL  JUN  84  OSURF-715107-84-4 
UNCLASSIFIED  AFOSR-TR-85-8092  AFOSR-83-0055  F/G  8713  NL 


microcopy  resolution  test 

national  bureau  Of  standards  i 


W  •  A  *  *  -  N  '  »  -  .  • 


AD- A 150  926 


IFOSR-TR. 


-009  2 


RF  Project  763420/715107 
Annual  Report 


RESPONSE  OF  SATURATED  SOILS 
TO  DYNAMIC  LOADING 


Ranbir  S.  Sandhu,  S.  J.  Hong  and  Baher  L.  Aboustit 
Department  of  Civil  Engineering 


DEPARTMENT  OF  THE  AIR  FORCE 
Air  Force  Office  of  Scientific  Research 
Bolling  Air  Force  Base,  D.C.  20332 

Grant  No.  AFQSR-83-0055 


DTIC 

ELECTE 
MAR  4  1986 


I 


85 


June,  1984  _ _ _ ^r-r 

ppiyfywimow  STATEMENT^ 

1  Apptc»*d  k*  P^tte 

\  Distribution  Pttlimited 


The  Ohio  State  University 
Research  Foundation 

1314  Kinnear  Road 
Columbus,  Ohio  43212 


oproved  for  public  roliati) 
3  036  ibution  unlimited. 


I.  REPOFIT  SECURITY  CLASSIf  (CATION 

Unclassified _ _ 


2»  SECURITY  CLASSIF  ICATION  AUTHORI T  V 


2t>  OECLASSIFICATlON/OOWNGRADING  schedule 


REPORT  DOCUMENTATION  PAGE 


1b.  RESTRICTIVE  MARKINGS 


3.  0ISTRI0UTION/AVAI  LABILITY  OF  REPORT 

Approved  for  Public  Release 
Distribution  Unlimited. 


4  PERFORMING  organization  REPORT  NUMBER  IE  I 

0SURF-715/0f-34-4 


S.  MONITORING  ORGANIZATION  REPORT  NUMBERISI 

AFOSR-TR.  85  -0  09  2 


6«  NAME  OF  PERFORMING  ORGANIZATION  Bb.  OFFICE  SYMBOL  7a.  NAME  OF  MONITORING  ORGANIZATION 

The  Ohio  State  Research  j  ",wUemkt*>  <- 

1  -"T 


laillliTir.lilhiii 


6c  AOORESS  (Cily.  Si  at*  and  /IP  Coda) 

Department  of  Civil  Engineering 
1314  Kinnear  Road 
Columbus,  OH  43212 


7b.  AOORESS  /City,  Stott  and  ZIP  Code > 


Ba.  NAME  OF  FUNOING/SPONSORING  Bb.  OFFICE  SYMBOL  B.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

ORGANIZATION  Air  Force  Of  appUeaUal 

Office  of  Scientific  Research  AFOSR/NA  AF0SR-83-0055 


8c  AOORESS  tCity.  State  and  ZIP  Code) 

Bolling  AFB,  DC  20332-6448 


10.  SOURCE  OF  FUNDING  NOS. 


PROGRAM 

element  no 


PROJECT 

NO 


11  TITLE  ilnclude  Security  Claeeitication)  KeS punse  UT  |  61102F 

Saturated  Soils  to  Dynamic  Loading  (Unclassified) 


WORK  UNIT 
NO 


13  PERSONAL  AUTHOR(S) 

Ranbir  S.  Sandhu  Soon  J.  Hon 


13a  TYPE  OF  REPORT  |  13b.  TIME  COVERED 

Annual 


16  SUPPLEMENTARY  NOTATION 


Baher  L.  Aboustit 


14  DATE  OF  REPORT  fVr  .  Mo..  Day t  15  PAGE  COUNT 

June,  1984  65 


ElSnaifil 


cosat i cooes 


18.  SUBJBCT  TERMS  (Continue  on  reverie  i f  neceuery  and  identify  by  block  number) 

Finite  Element  Methods 
Dynamics  of  Saturated  Soils 
Wave  Propagation 


19  ABSTRACT  -Continue  on  reverie  i f  neceuery  and  identify  by  block  numben 


The  transient  response  of  saturated  porous  soils  to  time  dependent  boundary  conditions 
is  analyzed.  Galerkin  finite  element  method  is  used  to  set  up  the  spatial  discretization 
of  E lot’s  equations  of  wave  propagation  through  linearly  elastic  fluid-saturated  porous 
medium.  Wilson's  B-y-5  algorithm  is  used  to  integrate  the  equations  of  motion.  The  procedur 
is  applied  to  several  one-dimensional  steady  state  and  transient  problems.  Excellent 
agreement  with  the  analytic  solutions  was  obtained  with  'proper'  selection  of  time- 
integration  parameters. 


30  O'STRIBUTION'AVAU.  ability  of  ABSTRACT 

unci  ASS1F1EO  UNLIMITED  ^  SAME  AS  RPT  Q  OTIC  USERS  O 


33*  name  of  RESPONSIBLE  INDIVIDUAL 


„ _ Lt.  Col.  Lawrence  D.  Hokanson 


DO  FORM  1473,  83  APR  e 


31  ABSTRACT  SECURITY  CLASSIFICATION 

Unclassified 


22b  TELEPHONE  NUMBER 
(Include  4r*«  Codei 

|22c  OFFICE  SvMBO 

202/767-4935 

'  AFOSR/NA 

EDITION  Op  *  JAN  73  IS 


ssru.*  *  v 


M  *■  ^  I  C  •  F  .\ 


T  p  Aok 


RESPONSE  OF 
SATURATED  SOILS 
TO  DYNAMIC  LOADING 


OSURF-7 151 07-84-4 


By 


Ranbir  S.  Sandhu 
S.J.  Hong 
Baber  L,  Aboustit 

Department  of  Civil  Engineering 


June  1984 


11*  TOBCS  OTFIC*  Of  SCI'EHT*** 

wrncs  0?  TKMIS’C.TTAL  TO  ^  _  ^  ^ t  %  „d  ls 

This  technic 

approved 

Dlstrlbutid  .-UwHed. 

M1TTKE*  J  »  ,  / 

Sl.t,  I..tal=al  IntorartloB  »»»**« 


AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEARCH 
Grant:  AF0SR-83-0055 


Geotechnical  Engineering  Report  No. 13 


The  Ohio  State  University  Research  Foundation 
1314  Kinnear  Road,  Columbus,  Ohio  43212 


FOREWORD 


The  investigation  reported  herein  is  part  of  the  research  project  at 
The  Ohio  State  University,  Columbus,  Ohio  supported  by  the  Air  Force  Of¬ 
fice  of  Scientific  Research  Grant  83-0055.  Lt.  Col.  John  J.  Allen  was 
the  Program  Manager  at  the  commencement  of  the  project.  Lt.  Col.  Law¬ 
rence  D.  Hokanson  was  the  Program  Manager  from  July  1,  1983  onwards. 
The  research  project  was  started  Feburary  1,  1983  and  is  continuing. 
The  present  report  documents  part  of  the  work  done  up  to  January  31, 
1984.  At  The  Ohio  State  University,  the  project  is  supervised  by  Dr. 
Ranbir  S.  Sandhu,  Professor,  Department  of  Civil  Engineering.  The  com¬ 
puter  program  modification  were  carried  out  by  Dr.  Baher  L.  Aboustit, 
Post-doctoral  Reseach  Associate.  The  analyses  reported  herein  were 
started  by  Dr.  Aboustit  and  completed  by  Mr.  Soon- jo  Hong,  graduate  stu¬ 
dent,  who  also  prepared  the  documentation  on  the  Ohio  State  University 
Computer  System.  The  Instruction  and  Research  Computer  Center  at  The 
Ohio  State  University  provided  the  computational  facilities. 


ABSTRACT 


The  transient  response  of  saturated  porous  soils  to  time  dependent 
boundary  conditions  is  analyzed.  Galerkin  finite  element  method  is  used 
to  set  up  the  spatial  discretization  of  Biot's  equations  of  wave  propa¬ 
gation  through  linearly  elastic  fluid-saturated  porous  medium.  Wilson's 
/3-y-9  algorithm  is  used  to  integrate  the  equations  of  motion.  The  proce¬ 
dure  is  applied  to  several  one-dimensional  steady  state  and  transient 
problems.  Excellent  agreement  with  the  analytic  solution  was  obtained 
with  'proper'  selection  of  time- integration  parameters.  There  is  appar¬ 
ent  need  for  developing  reliable  time-integration  procedures. 
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SECTION  I 


INTRODUCTION 

Biot  (2,3)  developed  the  field  equations  for  propagations  of  waves 
through  fluid-saturated  elastic  porous  media.  Lagrange  equation  of  mo¬ 
tion  was  set  up  assuming  the  existence  of  a  dissipation  function  in 
terms  of  the  relative  velocity  of  the  two  constituents.  By  assuming  the 
existence  of  a  strain  energy  density  function  for  the  mixture,  the  equa¬ 
tions  of  motion  were  written  in  terms  of  the  displacements  of  the  solid 
and  the  fluid.  Hsieh  and  Yew  (13)  combined  the  theory  of  mixtures  and 
Biot's  assumption  that  the  fluid  acts  over  the  pore  space,  to  provide 
similar  equations  to  Biot's  equation  (3)  but  with  variable  porosity. 
Garg  (10)  extended  the  theory  to  elastic-plastic  soils.  Ishihara  (15) 
examined  Biot's  work  in  terms  of  physical  constants  which  are  related  to 
the  compressibilities  of  individual  constituent  materials.  Prevost  (17) 
assumed  Newtonian  fluid  behavior  and  the  soil  to  be  elasto-plastic. 

Most  analytical  solutions  for  the  initial  boundary  value  problem  of 
wave  propagation  through  saturated  porous  media  provide  only  the  harmon¬ 
ic  component  of  the  solution.  Deresiewicz  (6,7)  obtained  solution  for 
the  reflection  of  plane  waves  (2)  and  Love  waves  at  a  free  plane.  Zienk- 
iewicz  et  al . ( 24 )  obtained  the  solution  for  a  half  space  subjected  to 
harmonic  tractions.  Garg  et  al.(ll)  presented  the  exact  transient  as 
well  as  the  steady  state  solutions  for  compressional  wave  propagation 
through  a  porous  elastic  solid  and  developed  finite  difference  proce- 


dures  for  numerical  solution.  Chakroborty  and  Dey  (4)  obtained  the  so¬ 
lution  for  Love  waves  in  saturated  media  underlain  by  heterogeneous 
elastic  media. 

Although  the  finite  element  method  has  been  extensively  used  for 
analysis  of  quasi-static  consolidation  problems,  e.g.  Sandhu  (21),  its 
application  to  the  dynamic  response  of  saturated  soils  is  still  in  its 
infancy.  Ghaboussi  and  Wilson  (12)  used  Sandhu  and  Pister's  (19)  ap¬ 
proach  to  construct  a  variational  principle,  of  the  Gurtin  type,  equiva¬ 
lent  to  Biot's  (3)  field  equations  including  initial  as  well  as  boundary 
conditions.  This  variational  principle  was  the  basis  for  a  finite  ele¬ 
ment  discretization,  in  which  the  Gurtin  type  approach  was  replaced  by 
the /8-  Y  -  9  method  for  the  time-domain  integration.  In  the  spatial 
discretization,  a  one-dimensional  element  was  used  with  nodal  values  of 
the  solid  displacement  and  the  relative  displacement  of  the  fluid  as 
generalized  coordinates.  The  problem  of  a  soil  layer  subjected  to  unit 
step  loading  was  analyzed.  The  problem  was  solved  by  two  methods  based 
on  lumped  mass  and  on  consistent  mass,  respectively.  Waves  appeared  to 
propagate  faster  when  consistent  mass  formulation  was  used.  This  was 
ascribed  to  inertial  coupling  in  this  formulation.  Ghaboussi  and  Wilson 
(12)  did  not  present  any  comparison  with  exact  solution.  Garg  (11)  con¬ 
sidered  this  unsatisfactory  as  the  results  from  the  two  approaches  dif¬ 
fered  consic’.,  ably.  Prevost  (17)  used  the  Galerkin  approach  for  spatial 
discretization  and  Hughes'  (14)  implicit-explicit  algorithm  for  time  do¬ 
main  integration.  Bilinear  four  noded  isoparametric  elements  were  used 
with  nodal  point  values  of  displacement  of  the  solid  and  the  fluid  as 
the  unknown  parameters.  Assumption  of  Newtonian  fluid  behavior  lead  to 


nonsymmetric  damping  matrices.  Two-dimensional  analysis  was  actually 
used  to  solve  a  one-dimensional  problem  with  finite  length.  For  rela¬ 
tively  large  permeability,  Prevost  compared  his  results  with  the  exact 
solution  for  a  solid  with  no  pores.  This  was  incorrect  because  non-po- 
rous  solid  represents  a  strong  coupling,  while  a  material  with  high 
permeability  is  weakly  coupled  (11). 

The  purpose  of  this  investigation  was  to  develop  an  effective  finite 
element  computer  program  based  on  the  Biot's  theory  and  evaluate  it 
against  available  analytical  solutions.  A  Galerkin  approach  was  used  to 
set  up  the  semi-discrete  matrices  spatially  and  the|3 -y-  8  algorithm  is 
used  for  the  time  domain  integration.  Nodal  values  of  the  solid  dis¬ 
placement  and  relative  displacement  of  the  fluid  with  respect  to  the 
solid  skeleton  are  used  as  the  unknown  quantities.  A  computer  program 
was  written  to  handle  one-dimensional  as  well  as  plane-strain  problems. 
One-dimensional  linear  element  and  four  node  and  eight  node  isoparame¬ 
tric  two-dimensional  elements  were  used.  The  development  is  essentially 
similar  to  Ghaboussi  and  Wilson's  (12).  The  program  was  applied  to  nu¬ 
merical  solution  of  several  one-dimensional  problems  for  quasi -static  as 
well  as  dynamic  problems  for  which  exact  solutions  are  available.  The 
results  showed  excellent  agreement  between  the  approximate  and  the  exact 
solution. 

SECTION  II  summarizes  the  governing  field  equations  following  Biot 
(3).  The  finite  element  development  is  given  in  section  III.  SECTION  IV 
describes  the  example  problems.  Results  of  the  numerical  solution  proce¬ 
dures  are  evaluated  in  SECTION  V. 


SECTION  II 


FIELD  EQUATIONS  GOVERNING  DYNAMICS  OF  FLUID-SATURATED 

ELASTIC  SOILS 

Biot's  (3)  equations  of  motion  for  an  elastic  porous  medium  saturated 
with  a  compressible  fluid  may  be  written  in  standard  indicial  notation 
as; 


[EijklU|c,l  +  aM(auk,k  +  wk,k>Sij],j  +  pfi  =  p“i  +  ““  P2-Wi  U) 

CH(auk,k 4  «k,k>],i 4  -7  Vi 4  -7  Vi +  72  Vi  +7  *i  <21 

where  u.f  ,  f^,  Eijkl  denote  the  cartesian  components,  respectively, 
of  the  solid  displacement  vector,  the  relative  fluid  displacement  vector, 
the  body  force  vector  per  unit  mass  and  the  isothermal  elasticity 
tensor.  P  is  mass  density  of  the  saturated  soil  and  P2  that  of  water  per 
unit  bulk  volume,  f,  K,  a,  M  are,  respectively,  the  porosity,  the  perme¬ 
ability,  the  solid  compressibility  and  the  fluid  compressibil ity.  The 
superposed  dot  implies  a  time  derivative.  All  the  functions  are  defined 
over  the  cartesian  product  R  x  [0,oo)  where  R  is  the  spatial  region  of 
interest  and  [0,«o)  is  the  positive  interval  of  time.  With  these  field 
equations,  we  associate  the  following  boundary  conditions. 

ui(t)  =  Uj(t)  on  Sjf  (3) 
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(4) 


ci  =rtJ"j  ’  (E(jkiuk,i  t0,r6ij)nj  *  li  on  s2i 
7T(t)  =  M(«uk>k  +  wk>|<)  =7T(t)  Ofl  S3  (5) 

w^t)  =  w^t)  on  S4  (6) 

where  ,  Sg.  are  complementary  subsets  of  the  boundary  S  of  the 
spatial  region  of  interest  and  so  are  S3,  S4.  The  initial  conditions  for 
the  problem  are  given  by: 

u(x,  0)  =  uq(x) 
u(x,  0)  =  uq(x) 
w(x,  0)  =  wQ(x) 
w(x,  0)  =  w  (x) 


SECTION  III 


FINITE  ELEMENT  FORMULATION 
3.1  SPATIAL  DISCRETIZATION 

Spatial  discretization,  APPENDIX  A,  of  the  governing  equations  for 
the  two-field  formulation  leads  to  the  following  matrix  equations; 
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The  development  of  these  vectors  and  matrices  is  given  in  Appendix  A. 
These  equations  assume  no  inherent  damping  in  the  system  as  a  whole. 
The  only  damping  component  is  associated  with  relative  motion. 

Ghaboussi  and  Wilson  (12)  introduced  damping  as  a  linear  combination 
of  the  stiffness  and  the  mass  matrix  in  the  following  form; 

Css=  al(Mss  ”  f2  Mff>  +  a2*Kss  "  0,2  Kff}  {8) 

where  a^  a2  are  constants  and  f,  of  have  been  defined  in  SECTION  II. 
The  structural  damping  matrix  is  a  linear  combination  of  the  mass  and 
the  effective  stiffness  of  the  soil.  Comparing  with  Equation  (1),  (K$s- 

°^Kff)u  corresponds  to  (Eij|(luk ,1 )  y  The  quantity  Mss  corresponds  to 
Pin  Equation  (1)  and  M^  to  ^/f2  in  Equation  (2).  Hence,  Mss-f^Mff 
corresponds  to  P1# 


3.2  TIME  DOMAIN  INTEGRATION 

The  discretized  equations  of  motion.  Equation  (7)  can  be  rewritten  as 

MU  +  CU  +  KU  =  R  (9) 

in  which 


Wilson's  /3  0  algorithm  was  used  to  integrate  the  equations  of  mot¬ 

ion.  In  summary,  the  algorithm  works  as  follows. 

The  displacements  and  velocities  at  time  tn+  0  j^t  can  be  expressed  in 

terms  of  U,  U,  il  at  time  t  as 

n 

Un+*  *  un  *  wt|j„  *  <1/2-®  «?<At)2  ii„  »/302(At)2  U^,  (11) 

'\t,  -  *  U->')  «At“„ 

i n  which  /J  and  V  .in*  Newmark  1 c  ucttu  i «*ot  •.  f  If*). 


.  t  I  I  u  I  lull 


1 1  1  •  mi  .t  i  l  *  if  >'  <  I  I  i  .in*l  M  ‘  t 


I .  I  *i 


I  •  I  •  i* .  t  ‘  *  J  /  I  •  >  * 


in  which 


*  1  y 

K  =  K  +  w  1  n  M  +  -  "  -  C 

Pm 


*  1  y 

^n+fl  =  ^r>4.a  +  "W  ~  H  Ma  +  '■■  Cb 
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Assuming  cubic  expansion  of  nodal  point  displacement  over  a  time 
interval  in  terms  of  displacement,  velocity  and  acceleration  at  time  t^, 
the  displacement,  velocity  and  acceleration  at  time  t^  +  At  can  be  ex¬ 
pressed  as 


W  >W  ♦  <l  -  pK  *  4t0„  +  l'2<>  -  f><  A‘)2 


^ntl  ■  TT - 5  (Un*.  •  un>  +  l1  -  ♦  <1  -  — >  tii„  (15) 

n  1  P^(M)2  n+®  n  002  n  2f}Q  n 


j"* 5  ife- (#»*.  ■  u-) -^r°-  * u  *5s> 5- 


-  *  -  *  _•>  *’  *•  •  '  ^  V  ..  •  *  .W*  .  k 


Various  values  for  b  and  0  can  lead  to  various  integration  schemes  (12). 
For  0=1,  the  stability  of  the/3-  y  method  has  been  examined.  Zienk- 
iewicz  (23)  gives  the  following  criteria  for  the  undamped  response  of  a 
linear  system. 

V  >  1/2 

0  >  l/4(l/2  +  Y)2  (16) 

1/2  +  0  +  Y>  0 


SECTION  IV 


EXAMPLE  PROBLEMS 

4.1  INTRODUCTION 

Several  examples  with  known  analytical  solution  were  solved  to  check 
the  validity  of  the  method  and  the  correctness  of  code.  These  examples 
included 

a.  Quasi-static  soil  consolidation 

b.  Dynamic  response  of  an  elastic  layer  of  single  material 

c.  Response  of  a  fluid-saturated  soil  layer 

All  the  problems  deal  with  compressional  wave  propagation  in  an  initial¬ 
ly  undisturbed,  homogeneous,  isotropic,  elastic,  porous  or  non-porous 
system  to  dynamic  loading.  The  system  for  c  above  was  subjected  to  spa¬ 
tially  uniform  surface  traction  q(t)  as  shown  in  Figure  1.  Solution  for 
non-porous  media  was  accomplished  by  prescribing  relative  displacement 
between  soil  and  fluid  to  be  zero. 

4.2  QUASI-STATIC  PROBLEM  OF  SOIL  CONSOLIDATION 

The  quasi-static  consolidation  problem  was  solved  by  setting  density 
and  damping  equal  to  zero.  The  example  problem  selected  for  cmparison 
was  the  one  solved  earlier  by  Sandhu  (20).  Excellent  agreement  was  ob¬ 
served  between  the  present  solution  and  the  one  obtained  (20)  using  a 


Figure  1:  Soi 
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quasi-  static  consolidation  analysis  computer  program  based  on  the  8-4 
element. 

4.3  DYNAMIC  RESPONSE  OF  AN  ELASTIC  BAR  OF  SINGLE  MATERIAL 

An  elastic  layer  under  spatially  uniform  excitation  applied  to  its 
surface  can  be  regarded  as  an  elastic  bar,  constrained  in  its  lateral 
dimensions  and  subjected  to  excitation  at  the  free  end.  A  representative 
segment  of  the  layer,  equivalent  to  a  laterally  restrained  bar,  is  shown 
in  Figure  2.  The  characteristics  used  in  the  finite  element  model  were; 


Total  length 
Number  of  nodes 
Number  of  elements 
Length  of  each  element 
Modulus  of  elasticity 
Poisson's  ratio 
Density 
Wave  velocity 
Time  interval 
Number  of  time  steps 


L  =  500  mm 
=  51 
=  50 
=  10  mm 

E  =  20000  kg/mm2 
v  =  0 

P  =  0.0008  kg.m.sec2/mm^ 

CQ  =  [(2M  +  A  )/P  ]1/2  =  5000  mm/msec 
At  =  0.002  msec 
=  50 


where  n,  A  are  Lame's  constants.  The  following  two  different  conditions 
at  the  free  end  were  considered. 

a.  STEADY  STATE  RESPONSE 
A  sinusoidal  loading  in  the  form 


q(*)  =  qQ  sin(407Tt) 


(19) 


was  applied  to  one  end  of  the  bar  with  the  other  end  fixed. 


b.  RESPONSE  TO  STEP  LOAD 

A  load  was  suddenly  applied  and  allowed  to  stay,  i.e. 

q(t)  =  qQ  H(t )  =PCQ  H  ( t )  (20) 

where  H(t)  is  the  Heaviside  function.  It  should  be  noted  that  the  load¬ 
ing  in  Equation  (20)  will  induce  a  unit  particle  velocity  (1).  Several 
values  for  £,  y,  0  were  used  to  obtain  the  optimum  solution.  The  ana¬ 
lytical  solution  for  this  problem  which  contains  both  the  harmonic  and 
the  transient  solution  is  given  in  reference  (9). 


4.4  RESPONSE  OF  A  FLUID-SATURATED  SOIL  LATER  (Garg's  Problem) 


The  third  example  considered  was  that  of  a  fluid-saturated  porous 
soil  layer.  The  finite  element  model  had  the  same  characteristics  as 
used  by  Garg  (11),  i.e.. 


Total  length 
Number  of  nodes 
Number  of  elements 
Length  of  each  element 
Modulus  of  elasticity 
Poisson's  ratio 
Mixture  mass  density 
Fluid  mass  density 
Porosity 

Fluid  compressibility 


L  =  50  cm 
=  51 
=  50 
=  1  cm 

E  =  0.2319xl012  dyn/cm2 
v  =  0.171 
P  =  2.3612  gm/cm2 
@2=  0.18  gm/cm2 
f  =  0.18 

M  =  0. 102xl012  dyn/cm2 
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Lower  bound  of  wave 

velocity  Cq  =  354875  cm/sec 

Upper  bound  of  first  kind 

wave  velocity  C+  =  358193  cm/sec 

Upper  bound  of  second  kind 

wave  velocity  C_  =  127941  cm/sec 

The  velocities  CQ,  C+>  C_  were  defined  by  Biot  (2)  and  Garg  (11).  Using 
the  notation  of  this  report,  these  quantities  are  given  by  the  following 
expressions  (11). 

C2  =  (  A+  IV-  +  a2  M )/p 
C2  =  (  A+  2m  +  M(a-f)2)/  P1 
C2  =  Mf2/  A, 

C22  =  Mf(a-f)/  Px 
C2n  =  Mf(o-f)/  P2 

2C+  =  Ci  +  C2  ±  L(C2  -  C2)2  +  4C22  C221]1/2 

Here,  =  P-  P2  is  the  bulk  mass  density  of  the  soil.  These  velocities 
are  applicable  to  one-dimensional  compressive  wave  propagation.  A  unit 
particle  velocity  for  each  phase  was  imposed  at  the  free  surface,  i.e.. 


u(L,t)  =  H(t ) 

(21) 

w(L,t)  =  0 

(22) 

Equation  (22)  implies  strong  coupling  at  the  free  surface.  In  reality, 
this  may  not  be  true.  However,  for  the  purpose  of  comparison  with  Garg, 
the  same  assumptions  were  made.  As  reported  by  Garg  et  el.  (11),  as 


permeability  K  — >  0,  the  relative  motion  between  the  two  constituents 
vanishes  and  phase  velocity  C  .  This  is  termed  "strong  coupling". 
In  this  case  the  material  behaves  as  a  single  continuum  whose  properties 
are  combination  of  those  of  the  two  constituents.  On  the  other  hand,  as 
K  — >oo  ,  the  coupling  between  the  two  constituents  vanishes  and  C  -*  C^. 
.  This  extreme  is  termed  "weak  coupling".  The  boundary  condition  given 
by  Equation  (21)  can  be  replaced  by  the  traction  boundary  condition, 
while  the  boundary  conditions  expressed  by  Equation  (22)  can  be  replaced 
by  the  displacement  boundary  condition,  i.e. 


q(L,t)  =  P CH ( t ) 

(23) 

w(L,t)  =  0 

(24) 

where  q(l,t)  is  the  traction  applied  to  the  free  surface.  The  numerical 
values  for  the  permeability  and  the  time  steps  corresponding  to  strong 
and  weak  couplings  were; 

a.  Strong  Coupling  (Low  Permeability) 

K  =  0.148x10"^  cm^/gm  sec 
At  =  1  micro  sec 
Number  of  time  steps  =  50 

b.  Weak  Coupling  (High  Permeability) 

K  =  0.148x10"^  cm^/gm  sec 


At  =  2.4  micro  sec 


Number  of  tme  steps  =  50 


4.5  RESPONSE  OF  A  FLUID-SATURATED  SOIL  LAYER  (Ghaboussi 's  Problem) 


The  fourth  example  was  the  problem  solved  by  Ghaboussi  and  Wilson 
(12).  They  did  not  compare  the  numerical  solution  with  any  analytical 
solution.  The  problem  considered  was  a  fluid-saturated  half  space  sub¬ 
jected  to  a  unit  step-loading.  The  finite  element  model  had  the  follow¬ 
ing  properties. 


Total  length 
Number  of  nodes 
Number  of  elements 
Length  of  each  element 
Modulus  of  elasticity 
Poisson's  ratio 
Mixture  density 
Fluid  density 
Porosity 

Fluid  compressibi 1 ity 

Solid  compressibi lty 

Wave  velocity  for  the 
mixture  (no  relative 
motion ) 

Coefficient  of  permeability 
Time  interval 


L  =  50  cm 
=  51 
=  50 
=  1  cm 

E  =  0.2319xl012  dyn/cm2 
v  *  0.171 
P  =  3  gm/cm3 
P2  =  0.18  gm/cm3 
f  =  0.18 

M  =  90xl0n  dyn/cm2 
a  ~  1 

C0  =  L(2H  +  x+  «2M)/P  ]1/2 
=  1.755895xl06  cm/sec 
K  =  0.19xl0"S  cm/sec 
At=  1  micro  sec 


Generally,  the  solid  density  of  soil  falls  between  the  range  of  2.0-  2./ 

3  10  2 

(gm/cm  )  and  the  compressibility  of  pure  water  is  2.0x10  (dyn/cm  ). 

Thus,  it  should  be  noted  that  mixture  density  and  the  fluid  rompr  1 1>  i  - 


lity  given  above  are  far  from  real  soil  properties,  but  were  taken  to 
match  the  problem  solved  by  Ghaboussi  and  Wilson  (12).  Reference  (12)did 
not  give  any  values  for  k  and  f.  Assigned  value  of  k  is  for  the  mixture 
of  very  fine  sand  and  silt  and  f  is  for  very  dense  state. 


SECTION  V 


RESULTS  OF  THE  ANALYSES 

5.1  QUASI-STATIC  PROBLEM  OF  SOIL  CONSOLIDATION 

The  results  for  the  present  analysis  and  those  obtained  by  Sandhu 
(20)  coincided  as  expected.  Detailed  results  of  the  analyses  were  as  re¬ 
ported  in  reference  (22)  for  the  problem  of  one-dimensional  consolida¬ 
tion  showed  using  Ghaboussi 's  4-4  element. 


5.2  DYNAMIC  RESPONSE  OF  AN  ELASTIC  LAYER  OF  SINGLE  MATERIAL 

Figures  3(a)  through  3(f)  and  Table  1(a)  thru  1(c)  illustrate  the 
stress  response  of  an  elastic  bar,  fixed  at  one  end  and  subjected  to  si¬ 
nusoidal  loading  at  the  free  end  at  time  stages  0.012,  0.038,  0.05, 

0.064,  0.088  and  0.10  seconds.  The  exact  solution  (9)  for  the  problem 
is,  for  applied  load  q(t)  =  qQsin(407rt), 

a.  Displacement: 

u(y,t)=  2  CG(6n)  -  G(*n+1)3  (25) 

n=0 

b.  Stress: 


a(y,t)=  -  -  £  CG'(0n)  +  G'(*n+1)] 

Cb  n=0 


(26) 
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a  NUMERICAL  SOLUTION 
®  EXACT  SOLUTION 


Stress  Distribution  in  an  Elastic  Layer 
Subjected  to  Sinusoidal  Load  at  Free 
Surface,  (a)  Times0.012  second. 


o  NUMERICAL  SOLUTION 
0  EXACT  SOLUTION 


Stress  Distribution  in  an  Elastic  Layer 
Subjected  to  Sinusoidal  Load  at  Free 
Surface,  (b)  Time=0.038  second. 


0  NUMERICAL  SOLUTION 
o  EXACT  SOLUTION 


Stress  Distribution  in  an  Elastic  Layer 
Subjected  to  Sinusoidal  Load  at  Free 
Surface,  (c)  Time=0.05  second. 


o  NUMERICAL  SOLUTION 
o  EXACT  SOLUTION 


Stress  Distribution  in  an  Elastic  Layer 
Subjected  to  Sinusoidal  Load  at  Free 
Surface,  (d)  Time*0.064  second. 


IlHglS] 


□  NUMERICAL  SOLUTION 
o  EXACT  SOLUTION 


Stress  Distribution  in  an  Elastic  Layer 
Subjected  to  Sinusoidal  Load  at  Free 
Surface,  (e)  Tirae=0.088  second 


0  NUMERICAL  SOLUTION 
o  EXACT  SOLUTION 


Figure  3 


Stress  Distribution  in  an  Elastic  Layer 
Subjected  to  Sinusoidal  Load  at  Free 
Surface.  ( f )  Time=0.100  second. 


TABLE  1 

Stress  Distribution  in  an  Elastic  Layer  of  Single  Material  Subjected  to 
Sinusoidal  Load  at  Free  Surface,  (a)  Time=  0.012  and  0.03  second 


t  =0 .012  sec 


FEM 

0.2349x10' 


0.6543x10 


0.1835x10 


0.4703x10 


0.107  xlO 


0.2058x10 


0.3073x10 


.3 


0.1231x10' 


0.3178 


0.9915 


EXACT 

0.0 


t=0.03 

sec 

FEM 

0.2033xl0'26 

0.0 

0.1104xl0"21 

0.0 

0.4833xl0'17 

0.0 

1 

0.1377x10"* 

0.0 

0.2132xl0"8 

0.0 

0.1297xl0"4 

0.0 

0.1581X10"1 

0.0 

0.7773 

0.7795 

0.8435 


0.368  -0.2454 


0.9823  -0.9423 


-0.2 


-0.9880 


a\ 


t=0.05 

sec 

t=0.Q64 

sec 

FEM 

EXACT 

I 

EXACT 

0.4820xl0-18 

0.0 

0.3508xl0-11 

0.0 

0.7783xl0-14 

0.0 

0.1262xl0"7 

0.0 

0.8188xl0“10 

0.0 

0.2138xl0-4 

0.0 

0.397  xlO"6 

_ 

0.0 

0.9445xl0-2 

0.0 

0.6025xl0-3 

0.0 

0.4798 

0.5878 

0.1343 

0.9235 

0.951 

0.9983 

0.9823 

0.2535x10-1 

0.0 

0.4748 

0.4828 

-0.9523 

-0.9510 

-0.6865 

-0.6845 

-0.5963 

-0.5878 

-0.911 

-0.9048 

0.5883 

0.5878 

-1.3225 


-1.2543 


0.997 


0.9980 
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TABLE  1 

Stress  Distribution  in  an  Elastic  Layer  of  Single  Material  Subjected  to 
Sinusoidal  Load  at  Free  Surface,  (c)  Time=  0.088  and  0.1  second 


a 

C+ 

II 

O 

« 

o 

CO 

00 

sec 

t=0.1 

sec 

FEM 

EXACT 

mmm 

EXACT 

0.4158xl0"5 

0.0 

0.2271 

0.1253 

0.6543xl0-2 

0.0 

1.0198 

0.9823 

0.8933 

0.9048 

0.5223 

0.4818 

0.689 

0.6845 

-0.6545 

-0.6845 

-0.4473 

-0.4818 

-0.918 

-0.9048 

-0.9913 

-0.9823 

0.1033 

0.1253 

-0.1448 

-0.1254 

0.9768 

0.9823 

0.893 

0.9048 

0.4865 

0.4818 

0.6815 

0.6845 

-0.6838 

-0.6845 

-0.4747 

-0.4818 

-0.901 

-0.9048 

0.99 


-0.994 


-0.998 


-0.1241 


-0.1254 


where 

G{t)=  Cb/E  Joq(t)  dt 

en  =  t  -  (2nL+L-y)/Cb 

<Pn  =  t  -  (2nL-L+y)/Cb 

dG 

G'({)  =  - 
dt 

G(0n)  =  0  where  0n<  0 

C2  =  E/P  (wave  velocity  of  bar) 

The  numerical  solution  was  based  on  values  of  qQ  =  -4,  in  addition  to 
those  given  in  section  4.3.  The  investigation  parameters  0,  y ,  0  were 

set  equal  to  0.25,  0.5,  1.0,  respectively.  These  values  satisfy  Equation 
(16).  In  Figures  3(a)  through  3(f),  we  see  good  coincidence  between  the 
finite  element  solution  and  the  analytical  solution  throughout  the  spa¬ 
tial  as  well  as  the  temporal  domains.  At  the  point  of  y/L=0.99,  the 
percentage  error  at  t=0.03  and  0.05  second  was  4.6  and  5.4 
.respectively.  But  at  the  other  time  stages  it  was  less  than  1.0*. 

Figures  4(a)  thru  4(c)  and  Table  2(a)  thru  2(c)  illustrates  the  dis¬ 
placement,  velocity  and  stresses  along  the  elastic  layer  under  unit  step 
load  at  the  free  surface,  at  time  t  =  0.08  sec,  i.e.  after  40  time  in¬ 
crements.  The  results  were  obtained  for  two  different  time  integration 


DISPLACEMENT 
.00  0.08 


y/L 

0 =0.167 

>3=0.25 

0.1 

0.4811x10-21 

4.3568x10-5 

0.2 

-0.1785x10-5 

6.051xl0-4 

0.3 

0.008636 

0.0091013 

0.4 

0.0194 

0.001909 

0.5 

0.02858 

0.02903 

0.6 

0.03945 

0.03905 

0.7 

0.04854 

0.04900 

0.8 

0.05948 

0.05903 

0.9 

0.06852 

0.06903 

1.0 

0.07949 

0.07896 

schemes  and  compared  with  analytical  solution  (1,9).  The  first  scheme 
used  9  =  1.0,  $  =  0.5,  Y  =  0.167,  and  the  second  0  =  1.0,  /3=  0.5,  Y- 
0.25.  It  is  shown  from  Figure  4  Table  2(a)  that  both  schemes  give  the 
almost  same  displacement  response  with  about  3  error  against  the  exact 
solution.  However,  )3  =  0.167  was  better  for  velocity  and  /9=  0.25  better 
for  stress  distribution.  Despite  this  fact,  large  error  around  wave 
front  was  still  observed  and  numerical  results  were  not  reliable.  The 
oscillatory  error  may  be  overcome  by  increasing  y  and  using  the  corre¬ 
sponding  value  of 0  from  Equation  (17)  as  we  shall  see  in  the  next  exam¬ 
ple.  For  this  problem,  one-dimensional  linear  element  (called  1-1  ele¬ 
ment)  was  used  and  CPU  time  on  the  AMDAHL  470/V8  Computer  was  15.18 
seconds  for  sinusoidal  loading  and  15.4  seconds  for  unit  step  loading, 
respectively. 

5.3  RESPONSE  OF  A  FLUID-SATURATED  SOIL  LAYER  (Garg's  Problem) 

In  order  to  investigate  the  effect  of  fluid-soil  interaction  on  wave 
propagation,  two  different  values  of  the  permeability  coefficient  were 
selected  to  approximate  "strong"  and  "weak"  coupling  extremes  described 
by  Garg  (11).  Figures  5(a)  and  5(b)  and  Table  3(a)  and  3(b)  show  the  ve¬ 
locity  for  both  the  solid  and  the  fluid  at  10cm  from  the  traction  bound¬ 
ary.  The  numerical  results  are  based  on  calculations  with  9  =  1,  /3  = 
0.6,  Y-  0.3025.  Reasonable  agreement  is  seen  between  the  finite  element 
and  the  analytical  solutions.  But  while  the  exact  solution  has  the  sharp 
disci ntinuity  in  the  wave  front,  numerical  was  diffused.  A  single  wave 
front  exists.  The  wave  is  propagating  with  velocity  =  C0,  and  the  solid 
velocity  is  the  same  as  the  fluid  velocity.  This  is  because,  for  this 
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a  NUMERICAL  SOLUTION 
o  EXACT  SOLUTION 


Figure  5:  (a)  Garg's  Problem  for  Strong  Coupling 

Velocity  History  of  the  solid  at  10cm 
from  the  Tractron  Boundary. 
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a  NUMERICAL  SOLUTION 
o  EXACT  SOLUTION 


Figure  5:  (b)  Garg's  Problem  for  Strong  Coupling 

Velocity  History  of  the  Fluid  at  10cm 
from  the  Traction  Boundary. 


TABLE  3 


problem,  relative  velocity  w  approaches  zero  and  the  two  constituents 
effectively  act  as  a  single  continuum. 

Figures  6(a)  and  6(b)  demonstrate  the  results  of  other  extreme  viz. 
"weak"  coupling.  The  values  of  0,  Y,  (i  were  the  same  as  for  the  "strong" 
coupling.  The  results  for  a  station  10  cm  away  from  traction  boundary  (y 
=  40cm)  are  quite  close  to  Gang's  analytical  solution  (11).  Existence 
of  two  wave  fronts  travelling  with  speeds  C_  and  C+  is  noticed. 

Pore  pressure  distribution  at  -nfferent  times  for  the  two  extremes  of 
"strong"  and  "weak"  coupling  are  plotted  in  Figure  7(a)  and  7(b).  A  sin¬ 
gle  phase  description  is  seen  in  Figure  7(a),  in  which  the  pressure  wave 
is  propagating  with  speed  (t) .  Figure  7(b)  clearly  demonstrates  the  exis¬ 
tence  of  two  waves  travelling  with  speed  C_  and  C^,  in  the  fluid  and  the 
solid,  respectively.  The  above  results  were  obtained  by  1-1  element. 
CPU  time  was  15.68  seconds  for  weak  coupling  analysis  and  15.4  seconds 
for  strong  coupling. 

5.4  RESPONSE  OF  A  FLUID-SATURATED  SOIL  LAYER  (GHABOUSSI * S  PROBLEM) 

Figure  8  presents  the  response  of  a  saturated  layer  to  a  unit  step 
loading.  Solution  of  this  problem  was  attempted  by  Ghaboussi  and  Wilson 
(12).  However,  their  results,  were  criticized  by  Garg  (11),  and  do  not 
agree  with  the  plot  in  Figure  8  (Table  5).  In  this  figure,  the  non-di¬ 
mensional  pore  pressure  7r/qQis  plotted  against  y*  =  y/(  KCQ).  For  a  =  1 
and  M  -Goethe  pore  pressure  should  be  equal  to  the  applied  traction. 
Ghaboussi  and  Wilson  (12)  reported  figures  which  do  not  match  Figure  8. 
With  1-1  element,  CPU  time  for  the  problem  was  12.35  seconds. 


GflRG  PROBLEM fWEAK  COUPLING) 


TIME  (MICRO  SEC) 


a  NUMERICAL  SOLUTION 
©  EXACT  SOLUTION 


Figure  6:  (b)  Garg's  Problem  for  Weak  Coupling; 

Velocity  History  of  the  Fluid  at  10cm 
from  the  Traction  Boundary. 
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TABLE  4 

Velocity  History  of  Garg's  Problem  at  10cm  from  the  fraction  Boundary 

for  Weak  Coupling. 


Time 

Solid  (cm/sec) 

Fluid  (cm/sec) 

(sec) 

FEM 

EXACT 

FEM 

EXACT 

0.0024 

-0.2009xl0"8 

0.0 

-0.5264xl0'7 

0.0 

0.0096 

0.1844xl0"6 

0.0 

-0.1150x10“^ 

0.0 

0.0192 

0.4323xl0-2 

0.0 

0.1890xl0"2 

0.0 

0.0312 

0.9575xl0“2 

1.0 

0.4078 

0.5 

0.0408 

1.0527 

1.0 

0.4428 

0.5 

0.0504 

1.0760 

1.0 

0.4834 

0.5 

0.060 

1.0856 

1.0 

0.4880 

0.5 

0.696 

1.0807 

1.0 

0.4720 

0.5 

0.0792 

1.0632 

1.0 

0.8662 

1.0 

0.0912 

1.0577 

1.0 

1.0260 

1.0 

0.1008 

1.0568 

1.0 

1.0475 

1.0 

0.1104 

1.0560 

1.0 

1.0555 

1.0 

0.120 

1.0562 

1.0 

1.0549 

1.0 

GflRG  PROBLEM (STRONG  COUPLING) 


0  TIME*48  MICRO  SEC 


A  1IME-72  MICRO  SEC 


+  7 1 ME  =  96  MICRO  SEC 


X  T I  ME* 1 20  MICRO  SEC 


Pore  Pressure  Distribution  at  Various 
Times  for  Garg's  Problem,  (a)  Strong 
Coupling. 


TABLE  5 


Pore  Pressure  (tt/q)  Distribution  of  Garg's  Problem 

0 


y/L 

Strong  Coupling 
(t=0.05  milli  sec) 

Weak  Coupling 
(t=0.12  milli  sec) 

0.01 

-0.1332x10"** 

0.5204xl0"2 

0.11 

-0.1296xl0-11 

0.8856X10"1 

0.21 

-0.9494xl0"9 

0.2026 

0.31 

-0.1544x10"® 

0.2031 

0.41 

-0.3558xl0"4 

0.2027 

0.51 

-0.2315x10"^ 

0.2027 

0.61 

-0.4647xl0_1 

0.2070 

0.71 

0.2345 

0.2747 

0.81 

0.2359 

0.2928 

0.91 

0.2344 

0.2928 

0.99 

0.2347 

0.2928 

TABLE  6 


SECTION  VI 


CONCLUSIONS 

Galerkin  finite  element  method  was  used  along  with  the  /9-y-8  algor¬ 
ithm  to  set  up  numerical  scheme  for  investigating  wave  propagation 
through  saturated  porous  media.  Several  problems  with  known  analytical 
solutions  were  analyzed.  Results  of  the  analyses  indicate  the  following: 

1.  The  integration  parameters  0,  yand  0  should  be  carefully  selected 
to  avoid  oscillatory  error. 

2.  The  scheme,  with  proper  selection  of  /3,  yand  0,  showed  exellent 
agreement  with  the  analytical  solutions. 

3.  The  numerical  and  analytical  results  show  the  importance  of  the 
role  of  permeability  in  single  or  double  phase  description  for  fluid- 
saturated  porous  media.  For  low  permeability,  there  is  little  relative 
motion  and  the  strong  coupling  on  single  material  description  would  be 
valid.  For  high  permeability,  the  two  phase  description  is  necessary. 

4.  The  computer  code  was  checked  only  for  one-dimensional  problems. 
Its  effectivenes  for  two-dimensional  (plane  strain)  problems  is  yet  to 
be  established. 

5.  Further  investigation  in  the  choice  of  damping  matrices  is  re¬ 
quired.  The  solution  process  for  a  few  idealized  problems  has  been 
checked,  but  the  assumptions  regarding  various  couplings  may  not  repre¬ 
sent  actual  soil  behavior. 


6.  The  computer  code  implemented  Ghaboussi  and  Wilson's  version  of 
Biot's  theory.  The  entire  fluid  mass  is  expected  to  be  in  relative  mo¬ 
tion.  In  other  theories,  an  interaction  mass  is  introduced.  This  would 
imply  a  "partial"  coupling  somewhere  between  the  "strong"  and  "weak" 
coupling  defined  by  Garg  (11).  Work  is  needed  to  quantify  the  coupling 
and  to  allow  for  this  in  the  computer  code  by  suitably  defining  mass  ma¬ 
trices. 

7.  The  computer  code  needs  to  be  extended  to  propagation  of  shear 
waves  and  Love  and  Rayleigh  waves.  Studies  are  needed  to  allow  for  re¬ 
flection  and  refraction  of  waves  at  interfaces  or  boundaries.  Dynamics 
of  nonhomogeneous,  anisotropic  and  nonlinear  soils  needs  to  be  investi- 
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Appendix  A 

SPATIAL  DISCRETIZATION 

In  this  appendix,  we  present  the  spatial  discretization  of  the  equa¬ 
tion  of  wave  propagation  through  saturated  porous  media  using  Galerkin's 
method. 
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yc e uk Tuk,i + « M6ijuk,k +  aM5ijwk,k]^nj ds+ 

S  v 

*/[Eijkluk,l  ^  *  a?  M  Vk.k  *  °  H  8ijwk,k 

V 

+  (P  5j  +  -pP^w. )  ^  dv  (A‘ 


where  are  the  components  of  the  outward  unit  normal  to  S.  The  total 
stress  tensor  is  defined  as  (4) 


C,,  =  E .  u.  ,  +  a  M  8,  -u.  L  +  o  M  84  .w 


'ij  i jkl  k,l 


ijuk,k  “  w1j"k,k 


(A— 4) 


Hence,  Equation  (A-3)  can  be  rewritten  as 


/V  i  dS  *  f’i  dv 

s  v 

■  /CE fjkluk,l  *H,i  *  «2  M  Vk.k  *1)  ♦ 0  M  sijwk,k 


+  (P^  +  J  p?vi. )  tfrM]  dv 


(A-5) 


From  the  boundary  conditions.  Equations  (3)  and  (4),  we  have 


\p  =  0  on  S, 


S21 


(A-7] 


Substituting  Equations  (A-6)  and  (A-7)  into  Equation  (A-5), 


r 


5 


■ 

pa 

i  ” 

I 

U 


i 


[  t,  0M  dS  ♦  /fi  0M  dv  =  /[Eijkluk#1  /j  ♦  (a2  M  5ijuk  k  0^ 

2i  v  v 

+  a  M  6ijWk^(c  0Mj  +  (P  +  -|-P2  w. )  0M]  dv  (A-8) 


Similarly,  a  weak  form  for  Equation  (2)  is 

0  ./[«(«  Uk>k  *  wkjk)_t  ♦  -ip2ft-  !<y i,  -  f4o2»(  -  d»  (A-9) 

V 

M 

where  0  is  a  test  function.  Equation  (A-9)  can  be  rewritten  as 
/[<M(a  uk>k  *  wk_k)  /)_,  -  M(a  ukjk  *  wkjk)  /, 

V 

*  rV i  -  7p25t  -  /?¥,•  -  ctjV*M  dv  ■ 0  <A-10) 

Noting  that  the  pore  pressure  is  given  by  (3)  as 


7T=  M(a  u.  +  w.  .)  ( A— 1 1 ) 

1  >  1  *  9  1 

and  using  Gauss  theorem  on  the  first  term.  Equation  (A-10)  can  be 
written  as 


yVni  <t>H  dS  jP2 f1  /  dv 
S  v  f 

’  >(a  uk  k  ♦  wk>„)  ♦  (±<*5,  ♦  $  8,  ♦  C^j)  ♦"]  d,  < A-12 ) 

V 

where  n^  are  components  of  the  outward  unit  normal  to  S.  Equations  (5) 
and  (6)  for  7r  •  0  are 
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0M  =  O 


(A-13) 


TTn  .  =  7Tn  . 
J  J 


(A-14) 


Hence,  Equation  (A-12)  may  be  written  as 


f  tfn.  0M  dS  ♦  f  j-  P2 f.  0M  dv 


=  /[M(a  uk>k  ♦  wkjk)  ♦  (j-  P2iii  +  ^  +  C^wj)  0M]  dv  ( A- 15 ) 


A. 2  GALERKIN  APPROXIMATION 

For  a  typical  element,  the  local  solid  and  relative  fluid  displace¬ 
ments  are  approximated  by 


■Ax)  uj 


(A— 16 ) 


0N(x)  w? 


( A— 17) 


where  t/rN  and  0N  are  interpolating  functions  defined  over  the  finite 
element  and  u^,  w*|*  are  the  values  of  solid  and  fluid  displacements  at 
the  node  N. 

Substituting  from  Equation  (A-16),  (A-17)  into  Equation  (A-8)  gives 


/  t^M  dS 


frc  /N  N  ,M  .  2  M  c  t/N  ,M  N 


,N  ,M  ..N 
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a  M  6^  w»  i/r  J  +  (P<AN  ii_.  +  i_  0N  (AM]  dv  ( A- 18 ) 


Define 


dv  =  Fsi 


f  t,  ds  ♦  ft.  / 

2i  v 

JL  i jkl  ^,1  v,j  u  ij  .,},kJ  ssik 


/“  H  5 f  J  A  ■  Ksf1k 


dv  =  M1 


NM 

ss 


I 


f «,  /  /  dv  -  M™ 


Substituting  Equation  (A-19)  into  (A-18), 


kNM  N  NM  N  .  mNM  ..N  .  mNM  aN  CM 
Kssik  uk  +  Ksfik  wk  +  Mss  ui  +  Msf  wi  =  Fsi 


(A-19) 


(A— 20) 


Similarly,  substituting  from  Equation  (A-16), (A-17)  into  Equation  (A-15) 

f  0M7rn.  dS  +y  y  P2f.  9>M  dv 
S3  v 

■/"«« 

V  V  V 

+  /  ^  P2^N{ii  dv  *f  cij  ^N*jN  dv  (A-21) 


Define 


if  1 


